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ABSTRACT 

We have measured hne-of-sight velocity profiles (vPs) in the EO galaxy NGC 6703 
out to 2.6i?e- Comparing with the VPs predicted from spherical distribution functions 
(dps), we constrain the mass distribution and the anisotropy of the stellar orbits in 
this galaxy. 

We have developed a non-parametric technique to determine the df f{E, L"^) 
directly from the kinematic data. We test this technique on Monte Carlo simulated 
data with the spatial extent, sampling, and error bars of the NGC 6703 data. We find 
that smooth underlying DFs can be recovered to an RMS accuracy of 12% inside three 
times the radius of the last kinematic data point, and the anisotropy parameter j3{r) 
to an accuracy of 0.1, in a known potential. These uncertainties can be reduced with 
improved data. 

By comparing such best-estimate, regularized models in different potentials, 
we can derive constraints on the mass distribution and anisotropy. Tests show that 
with presently available data, an asymptotically constant halo circular velocity vq can 
be determined with an accuracy of ± < 50kms^^. This formal range often includes 
high-tiQ models with implausibly large gradients across the data boundary. However, 
even with extremely high quality data some uncertainty on the detailed shape of the 
underlying circular velocity curve remains. 

In the case of NGC 6703 we thus determine the true circular velocity at 2.6i?e 
to be 250 ± 40kms~^ at 95% confidence, corresponding to a total mass in NGC 6703 
inside 78" (13.5 /i^o^kpc, where /150 = ifo/50km/s/Mpc) of 1.6-2.6 x W^h^^ Mq. No 
model without dark matter will fit the data; however, a maximum stellar mass model 
in which the luminous component provides nearly all the mass in the centre does. In 
such a model, the total luminous mass inside 78" is 9 x 10^° Mq and the integrated 
B-band mass-to-light ratio out to this radius is = 5.3 — 10, corresponding to a 
rise from the center by at least a factor of 1.6. 

The anisotropy of the stellar distribution function in NGC 6703 changes from 
near-isotropic at the centre to slightly radially anisotropic (/3 — 0.3 — 0.4 at 30", 
(3 = 0.2 — 0.4 at 60") and is not well-constrained at the outer edge of the data, where 
/3 = — 0.5 — 1-0.4, depending on variations of the potential in the allowed range. 

Our results suggest that also elliptical galaxies begin to be dominated by dark 
matter at radii of ~ 10 kpc. 

Key words: stellar dynamics - dark matter - galaxies: elliptical and lenticular - 
galaxies: kinematics and dynamics - galaxies: individual - line profiles 
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1 INTRODUCTION 

Current cosmological models predict that, similar to spiral 
galaxies, elliptical galaxies should be surrounded by dark 
matter haloes. The observational evidence for dark matter 



© 0000 RAS 



2 O.E. Gerhard, G. Jeske, R.P. Saglia, R. Bender 



in ellipticals is still weak, however. In a few cases masses 
have been determined from X-ray observations (e.g., Awaki 
et al. 1994, Kim & Fabbiano 1995) or HI ring velocities 
(Franx, van Gorkom & de Zeeuw 1994). In others it has 
been possible to rule out constant M/L from extended ve- 
locity dispersion data (Saglia et al. 1993), from absorption 
hne profile measurements (CaroUo et al. 1995, Rix et al. 
1997), or from globular cluster or planetary nebula velocities 
(e.g., Grillmair et al. 1994, Arnaboldi et al. 1994). Gravita- 
tional lensing statistics (Maoz & Rix 1993) and individual 
image lens separations (Kochanek & Keeton 1997) favour 
models with extended dark matter haloes around ellipticals. 
Despite of this, the detailed radial mass distribution in ellip- 
tical galaxies remains largely unknown. Similarly, although 
we know from the tensor virial theorem that giant ellipti- 
cals are globally anisotropic (Binney 1978), their detailed 
anisotropy structure is only poorly known. 

The origin of this uncertainty is a fundamental de- 
generacy - in general, it is impossible to disentangle the 
anisotropy in the velocity distribution and the gravitational 
potential from velocity dispersion and rotation measure- 
ments alone (Binney & Mamon 1982, Dejonghe & Merritt 
1992). Tangential anisotropy, for example, can mimic the 
presence of dark matter. Recent dynamical studies have in- 
dicated, however, that the anisotropy of the stellar distribu- 
tion function (DF) is reflected in the shapes of the line-of- 
sight velocity profiles (VPs) in a way that depends on the 
gravitational potential (Gerhard 1993, G93; Merritt 1993). 
These papers argued that the extra constraints derived from 
the VP measurements may be enough to break the degener- 
acy and determine the mass distribution. 

If this is correct, it provides a new method to investi- 
gate the properties of the dark matter haloes around ellipti- 
cal galaxies at intermediate radii: VPs can now be estimated 
from high-quality absorption line measurements out to ~ 3 
effective radii. Dynamical models are then used to disentan- 
gle the effects of orbital anisotropy and potential gradient 
on the VP shapes. In this paper, we implement these ideas 
for the analysis of real data, analyzing the EO galaxy NGC 
6703. This study is part of an observational and theoretical 
program aimed at understanding the mass distribution and 
orbital structure in elliptical galaxies. Preliminary accounts 
of this work have been given in Jeske et al. (1996) and Saglia 
et al. (1997a). 

We have obtained long-slit spectroscopy for NGC 6703, 
and have measured VPs to ~ 2.6-Re with the method of Ben- 
der (1990). The results are quantified by a Gauss-Hermite 
decomposition (G93, van der Marel & Franx 1993) as de- 
scribed by Bender, Saglia & Gerhard (1994; BSG). These 
observations are described in Section 2. In Section 3 we use 
simple dynamical models to describe the variation of the 
VP shapes with anisotropy and potential, generalizing the 
results of G93 for scale-free models. These models, taken 
from a systematic study of the relation between DF and VPs 
in spherical potentials (Jeske 1995), are described in Ap- 
pendix A. In Section 4 we develop a non-parametric method 
for inferring the DF and potential from absorption line pro- 
file measurements. Tests on Monte Carlo generated data are 
used to determine the degree of confidence with which the 
DF and potential can be inferred from real data. In Section 
5 we analyse the kinematic data for NGC 6703. Compari- 
son with the dynamical models from Jeske (1995) already 



shows that no constant-M/L model will fit the data. The 
non-parametric method developed in Section 4 is then used 
to derive quantitative constraints on the mass distribution 
and anisotropy of this galaxy. Finally, Section 6 presents a 
discussion of the results and our conclusions. 



2 OBSERVATIONS OF NGC 6703 

NGC 6703 is an EO galaxy at a distance D = 36 Mpc 
(Faber et al. 1989) for Hq = 50 km/s/Mpc. From a B- 
band CCD frame taken at the prime focus of the 3.5m 
telescope on Calar Alto and kindly provided by U. Hopp 
we have measured the inner surface brightness profile us- 
ing the algorithm by Saglia et al. (1997b); this follows an 
7?^^* law with 7?e = 30" = 5.2 /iJq^ kpc, or a Jaffe model with 
rj — 46.5" — 8.1/iJo^kpc, with small residuals (Fig. 1). A 
3% increase of the sky value reduces the measured Jaffe ra- 
dius to 35.5", a 1% decrease of the sky values increases it 
to 54". Isophote shapes deviate little from circles (e<0.05, 
I 04/0 |< 0.005) and show smaU twisting {APA ~ 10°). 
From the Jaffe profile fit we derive a fiducial (calibrated and 
corrected for galactic absorption following Faber et al. 1989) 
Mi3 = -21.07, or luminosity Lis =4.16 x IQ^'^h^^ Lq^b- Note 
that the values of and Mb derived here are slightly larger 
than those {Re = 24", Mb = -20.79) given by Faber et al. 
(1989). 

The spectroscopic observations were carried out in Oc- 
tober 1994, May 1995 and August 1995 with the 3.5-m tele- 
scope on Calar Alto, Spain. In all of the runs the same setup 
was adopted. The Boiler & Chivens longslit twin spectro- 
graph was used with a 1200 line mm~^ grating giving 36 
A/mm dispersion. The detector was a Tektronix CCD with 
1024 X 1024 24^m pixels and the wavelength range 4760-5640 
A. The instrumental resolution obtained using a 3.6 arcsec 
wide slit was 85 km/s. We collected 1.5 hrs of observations 
along the major axis of the galaxy, 4 hrs of observations per- 
pendicular to the major axis and shifted to the North-East 
of the center by 24 arcsec (0.8i?e), and 13 hours of obser- 
vations perpendicular to the major axis and shifted to the 
North-East of the center by 36 arcsec (1.2i?e). Spectra taken 
parallel to the minor axis and shifted from the center allow 
at the same time a good sky subtraction and the symmetry 
check of the data points. 

The analysis of the data was carried out following the 
steps described by BSG. The logarithmic wavelength cali- 
bration was performed at a smaller step (Au = 30 km/s) 
than the actual pixel size (~ 50 km/s) to exploit the full ca- 
pabilities of the Fourier Correlation Quotient method. A sky 
subtraction better than 1% was achieved. The heliocentric 
velocity difference between the May 1995 and the Septem- 
ber 1994 - August 1995 frames was taken into account before 
coadding the observed spectra. The spectra were rebinned 
along the spatial direction to obtain a nearly constant sig- 
nal to noise ratio larger than 50 per resolution element. The 
effects of the continuum fitting and instrumental resolution 
were extensively tested by Monte Carlo simulations. The 
residual systematic effects on the values of the /13 and hi 
parameters are expected to be less than 0.01. The resulting 
fitted values for the folded velocity v , velocity dispersion 
a, /13 and /14 profiles are shown in Fig. |^ as a function of 
the distance from the center, reaching ~ 2.67?e. The v and 
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Figure 1. (a) Residuals of a Jaffe model fit to photometry for 
NGC 6703. Full line: in surface brightness; dotted line: in the 
curve of growth. The vertical dotted line marks the scaling Jaffe 
radius rj. (b) Folded mean velocity, (c) velocity dispersion, (d) 
?t3, and (e) ?t4 profiles. Crosses and filled circles refer to the two 
sides of the galaxy and the major axis spectrum. The small dots 
refer to the unrebinned spectrum (see text). Open and filled tri- 
angles refer to the two sides of the galaxy and the spectrum taken 
parallel to the minor axis and shifted 24 arcsec from the center. 
Open and filled squares refer to the two sides of the galaxy and the 
spectrum taken parallel to the minor axis and shifted 36 arcsec 
from the center. 



G profiles are folded antisymetrically with respect to the 
center for the major axis spectrum, symmetrically with re- 
spect to the major axis for the spectra parallel to the minor 
axis. Template mismatching was minimized by choosing the 
template star which gave the minimal symmetric profile 
derived along the major axis of the galaxy. The systematic 
effect due to the residual mismatching on the derived hi val- 
ues estimated from the remaining symmetric components is 
less than 0.01. 

The galaxy shows very little rotation (~ km/s for 
J? < 7?e, « 20 - 30 km/s for i? > iie). The (cylindrical) ro- 
tation measured parallel to the minor axis is slightly larger 
(~ 35 km/s) along the 24 arcsec shifted spectrum, but con- 
sistent with the peak velocity reached along the major axis 
at i? « 22 arcsec. This is shown by the velocities derived 
from the unbinned major axis spectra (dots in Fig. The 
velocity dispersion drops from the central « 190 km/s to 
~140 km/s at Re/2, slowly declining to about 110 km/s in 
the outer parts. The and /14 values are everywhere close 
to zero. The error bars are determined from Monte-Carlo 
simulations. Noise is added to template stars (rebinned to 
the original wavelength pixel size) broadened following the 
observed values of a and /14, matching the power spectrum 
noise to peak ratio of the galaxy spectra. The accuracy of 
the estimated error bars (the r.m.s. of 30 replica of the data 



points) is about 20 % (determined from the scatter of the 
estimated signal to noise ratios). 

Data points at the largest distances for the different 
data sets have lower signal to noise ratio than the mean 
and therefore have larger error bars. In addition, these data 
points are expected to suffer more from the systematic ef- 
fects due to the galaxy light contamination of the sky sub- 
traction(see discussion in Saglia et al. 1993). They are in any 
case consistent within the error bars with the more accurate 
values derived from the other available spectra. 

The observed scatter is sometimes slightly larger than 
expected from the error estimates. This excess could be real 
and due to the faint structures apparent in an unsharped 
mask image of the galaxy. In particular this applies to the 
asymmetries observed in the /14 profile in the central 5 arc- 
sec. The negative /14 values detected for the first data points 
of the 24 arcsec shifted spectrum are also real. They are de- 
tected in the unbinned major axis spectra at « 22 arcsec 
(dots in Fig. 



3 VELOCITY PROFILES IN SPHERICAL 
GALAXIES 

To better understand the relation between VP-shape, 
anisotropy, and gravitational potential, we have constructed 
a large number of anisotropic models for spherical galaxies 
in which the stars follow a Jaffe (1983) profile. The gravi- 
tational potential was taken to be either that of the stars 
(self-consistent case) , or one with everywhere constant rota- 
tion speed ('halo potential'). The latter case corresponds to 
a mass distribution with a dark halo which has equal density 
as the stars at r~0.4rj, and equal interior mass at r~rj, 
where r j is the scale radius of the Jaffe model. Anisotropic 
quasi-separable distribution functions (DFs) g{E)h{E,L^) 
were calculated by the method of Gerhard (1991; G91), but 
contrary to G93 the circularity function (which specifies the 
distribution of angular momenta on energy shells) was al- 
lowed to vary with energy. These models include DFs in 
which the anisotropy changes radially, from tangential to 
radial or vice-versa, from isotropic to radial to tangential, 
etc. (see Fig. |l^ ). For comparison, we have also constructed 
families of Osipkov (1979) - Merritt (1985) models which be- 
come strongly radially anisotropic beyond a certain radius. 
Details of this model data base and the properties of their 
VPs are given in Appendix A and in Jeske (1995); here we 
only give a brief summary relevant for the comparison with 
NGC 6703. 

The shapes of the observable VPs are most sensitive to 
the anisotropy of the DF, but depend also on the poten- 
tial (G93). For rapidly falling luminosity profiles, the VPs 
are dominated by the stars at the tangent point. Then ra- 
dially (tangentially) anisotropic DFs lead to more peaked 
(more flat-topped) VPs than in the isotropic case; in terms 
of the Gauss-Hermite parameter ^4, this corresponds to 
/i4>(/i4)iso and /i4<(ft4)iBo, respectively (Figs. 8,9 in G93). 

Fig. ^ shows that these trends are also seen in the 
present models in which both the luminosity density and the 
anisotropy change with radius. An increase in radial (tan- 
gential) anisotropy at intrinsic radius r is accompanied by 
an increase (decrease) of /14 at projected radius R~r. The 
correspondence is strongest in the models' outer parts, but 
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Figure 2. Fitted projected velocity dispersion tTfjt, anisotropy parameter /3, and VP-parameter h4 for representative Jaffe models in 
self-consistent (left) and halo potential (right). The models shown are radially and tangentially anisotropic models constructed with the 
method of G91 (dashed lines), the isotropic model (solid), and two Osipkov-Merritt models (dotted lines). Note that while /3 is a function 
of three-dimensional radius r and iTfit, /i4 are observed quantities depending on projected radius R, there is a close correspondence 
between features in these profiles. See text. 



is also seen to a lesser extent in the centre of a JafTe model 
where p{r) oc - contrary to a homogeneous core where 
radial orbits lead to broadened VPs (Dejonghe 1987). Quan- 
titatively the correspondence depends also on the anisotropy 
gradient. Osipkov-Merritt-models show a reversal of this 
trend near their anisotropy radius Va due to the large num- 
ber of high-energy radial orbits all turning around near r-^; 
this leads to flat-topped VPs in a small radius range near 
Ta. However, the properties of these models are extreme and 
they are in general not very useful for modelling observed 
VPs. 

Fig. ^ and Figs. 8,9 in G93 also show that as the mass 
of the model at large r is increased at constant anisotropy, 
both the projected dispersion and hi increase. Increasing 
P at constant potential, on the other hand, lowers a and 



increases /14. This suggests that by modelling a and /14 both 
mass M{r) and anisotropy /3(r) can in principle be found. 
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4 MODELLING ABSORPTION LINE PROFILE 
DATA 

Having seen the effect of anisotropy and potential variations 
on the line profile parameters, we now proceed to construct 
an algorithm by which the distribution fianction and poten- 
tial of a spherical galaxy can be constrained from its ob- 
served a and /14-profiles. Such absorption line profile data 
contain a subset of the information given by the projected 
distribution function N{R,V\\), which in the spherical case 
is related to the full DF by 



(1) 



Jdz J J dvxdvyfl^v, 



2.12. I,/ 1 2 2i 

-I- -^Vt + $(r),r Vt] 



Here the position on the sky is specified by i? = {x,y); 
R = \R\. Velocities in the sky plane are denoted by {vx,Vy), 
z and v\\ are the line-of-sight position and velocity, and Vr 
and Vt arc the intrinsic radial and tangential velocities. In 
spherical symmetry, the DF f{E,L^) is a function of energy 
and squared angular momentum only. 

Notice from eq. (1) that the projected DF depends lin- 
early on /, but non-linearly on the potential (l>(r). Thus 
/ will generally be easier to determine from iV(i?, uy) than 
Moreover, while considerations like those in the last sec- 
tion do suggest that, in spherical symmetry and for positive 
f{E,L^), both the DF and the potential can be determined 
from N{R,v\\), there is no theoretical proof that this is in 
fact true. We only know that in a fixed spherical potential 
the DF is uniquely determined from N{R,V\\) (Dejonghe & 
Merritt 1992). For these reasons we have found it useful to 
split our problem into two parts (see also Merritt 1993): 

(1) We fix the potential $, and from the photometric 
and kinematic data determine the "best" DF / for this po- 
tential. Because in practice the surface brightness (SB) pro- 
file is much better sampled than the kinematic observations 
and also has smaller errors, we treat it separately and de- 
termine the stellar luminosity density j{r) at the beginning. 
The kinematic data are then used to determine the "best" / 
for given j{r) and $(r), by approximately solving equation 
(1) as a linear integral equation. 

(2) We then vary $ to find that potential which allows 
the best fit overall. At present, it is not practical in step (2) 
to attempt to determine the potential non-parametrically. 
Rather, we choose a parametrized form for and find the 
region in parameter space for which the "best" DF as deter- 
mined in step (1) reproduces the data adequately. 

In view of the modelling of NGC 6703 in Section 5, we 
have considered the following family of potentials, includ- 
ing a luminous and a dark matter component: The stellar 
component is approximated as a JafFe (1983) sphere, with 
scale-radius rj and total mass Mj, so that 

^ , , GMj , r 

*L (r) = In . (2) 

rj r + rj 

The dark halo has an asymptotically flat rotation curve, 

Vc{r)=vo-^ J' , (3) 
so that its potential is 

^H{r) = ^v^ln{r' + rl). (4) 



This is specified by the asymptotic circular velocity vo and 

the core radius Vc. Both the luminous and dark halo compo- 
nents can be modified when needed, and need not be analytic 
functions. 

In testing our method below, we use parameters 
adapted to NGC 6703. This galaxy is well-fit by a Jaffe 
profile (Section 2), so rj is known. This leaves three free 
parameters, the mass Mj or mass-to-light ratio T of the 
stellar component, and the halo parameters rc and vo- If 
one assumes that the central kinematics is dominated by 
the luminous matter, T can be determined. Then only the 
two halo parameters Tc and vq are free. The assumption of 
maximum stellar T is similar to the maximum disk assump- 
tion in spiral galaxies. 

In any of the potentials specified by eqns. (2)-(4) we 
determine the DF by the algorithm described in Section 4.1 
below. To assess the significance of the results obtained, we 
test the algorithm on Monte Caxlo-generated pseudo data in 
Section 4.2. For kinematic data with the spatial extent and 
observational errors such as measured for NGC 6703, the 
algorithm recovers a smooth spherical DF ~ 70% of the time 
to an RMS level of ~ 12%, taken inside three times the radius 
of the outermost kinematic data point. In Section 4.3, we 
investigate the degree to which the gravitational potential 
can be constrained from similar data. 



4.1 Recovering / from a and /i4, given $ 

As discussed above, the projected distribution function 
N{r, V\\) suffices to determine DF f{E, L^) uniquely. In prac- 
tice, however, only incomplete and noisy data are available 
in place of A'^(r, V||), and contrary to the two-dimensional 
function N{r,V\{), the observed o{Ri) and h4,{Ri) contain 
only one-dimensional information. This suggests that we can 
hope to recover only the gross features of / from such kine- 
matic data. Indeed, the anisotropy parameter /3(r) seems 
to be essentially fixed from accurate measurements (e.g.. 
Figs. 8,9 in G93). Local fiuctuations in the DF will be inac- 
cessible, but as we will show, smooth DFs can be recovered 
with reasonable accuracy from presently available data. 

To solve the inversion problem, we first compute a set of 
self-consistent models fk{E,x) for the stellar density j{r), 
in the fixed potential "l>(r). The fk{E,x) are models of the 
kind discussed in Section 2 and Appendix A; E and x are 
the energy and circularity integrals of the motion. Then we 
write the DF as a sum over these "basis" functions: 



/= ^ akfk{E,x). 



(5) 



k=l,K 



We do not need to use a doubly infinite, complete set of 
basis functions because most of the high-frequency struc- 
ture represented by the higher-order basis functions in such 
a set will be swamped by noise in the observational data. 
It is sufficient to choose the number of basis functions, K, 
and the fk{E,x) themselves such that the data can bo fit- 
ted with a mean x'^ — 1 per data point. We have found it 
advantageous to use the isotropic model plus tangentially 
anisotropic basis models, because with these the anisotropy 
of the final composite DF (5) can be varied in a more local 
way than with radially anisotropic components. Since the 
ttfe can be negative, it is of course no problem to generate a 
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Figure 3. Anisotropy parameter /3 for a subset of basis functions 
/fc used in tlie self-consistent JafTe potential (top) and in a mixed 
Jaffe plus halo model (bottom). The top full line in each panel 
shows the isotropic model. All other basis models are tangentially 
anisotropic. Two values of asymptotic anisotropy as r — > oo are 
used (full and dashed lines). For illustration, the long-dashed lines 
show the /3-profiles of the best— fitting DFs derived with these 
bases from the NGC 6703 data, in both potentials. 



radially anisotropic DF from the isotropic model plus a set 
of tangential basis models. 

In Fig. 1^ we illustrate the basis models we have used, 
by plotting radial profiles of the anisotropy parameter /3(r) 
for a subset of them. For these basis functions /3{r) can be 
regarded as a measure of the extent of the DF in circularity 
X, at the energy E = $(r). The figure shows that the basis 
resolves three steps in x in the limits r — > and r —> oo; 
at intermediate energies, it has much finer resolution. The 
majority of the number of functions is spent on resolving 
the energy dependence, i.e., on placing the main gradient 
zones of the functions fk in the {E, x) plane on a relatively 
dense grid in energy. The main difference to using power law 
components (Fricke 1952, and, e.g., Dejonghe et al. 1996) is 
that our basis models are already reasonable in the sense 
that they are viable dynamical models for the density in 
question, and that the superposition is used only to match 
the kinematics. We have typically used K ^ 20 such func- 
tions; this proved sufficient even for analysing pseudo data 
of much better quality than we have for NGC 6703. 

Each of the fk{E,x) reproduces the stellar density dis- 
tribution j{r); so the ak satisfy 



a.k = 1- 



(6) 



Moreover, the coefficients ak must take values such that 

ak.fk{E,x) > (7) 

k = l,K 

everywhere in phase-space. In practice, these positivity con- 
straints are imposed on a grid in energy and circularity 



{E — Ei,x — Xj). Subject to these constraints the are 
to be determined such that the kinematics predicted from 
the DF (5) match the observed kinematics in a minimum 
X^-sense. 

The comparison between model and data is not entirely 
straightforward, however. Since the measured {v, a, /13, h^) 
are obtained by Biting to the line profile (the observational 
analogue of the projected DF), they unfortunately depend 
non-linearly on the galaxy's underlying DF. Thus we cannot 
use the observed v and cr in a linear least squares algorithm 
to determine the at from the data - they cannot be writ- 
ten as moments of /. In the fitting process, they have to be 
replaced by quantities that do depend linearly on the DF. 
Moreover, the error bars for these new quantities generally 
depend not only on the observed error bars of v and a, but 
also on the errors and the values of the line profile param- 
eters /13 and /i4. They must therefore be determined with 
some care. 

We have investigated several schemes along these lines. 
The following seemed to do best in recovering a known un- 
derlying spherical DF from pseudo data. From the measured 
{a, hi), we compute an approximation to the true velocity 
dispersion (second moment) a^{Ri), by integrating over the 
line profile (for negative h^, until it first becomes negative). 
We also evaluate a new set of even Gauss-Hermite moments 
s„{Ri;(j) from the data, using fixed, fiducial velocity scales 
cr{Ri) (G93, Appendix B). We have found it convenient to 
take for these a{Ri) the velocity dispersions (Tiso(-Ri) of the 
isotropic model with the galaxy's stellar density, in the cur- 
rent potential <I>(r). 

The velocity profile moments of the basis function mo- 
dels are transformed to the same velocity scales a{Ri). The 
moments a^{Ri) and Sn{Ri\a) of the composite DF then 
depend linearly on the corresponding moments of the basis 
function models. For a regularized model they are smooth 
functions of R, of which the (noisy) observational moments 
are assumed to be a random realisation within the respective 
errors. 

To determine the best-fitting coefficient ak we minimize 
the sum over data points i of all 



2 
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wl{Rr) 



■ w„(Ri) 



akal{Ri) 



,{R,;a) 
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k^l,K 



akS 



(fc) 



{R^ 



(8) 



(9) 



for n = 2,4. These equations make use of the fact that all 
the fk are self-consistent models for the same j(r), so that 
all surface density factors HkiRi) ~ f^{Ri) cancel. 

To determine the weights Wa and Wn, we have done 
Monte Carlo simulations studying the propagation of the 
observational errors Act and Ah^. Based on the results of 
these simulations, we have chosen 

w-\R,) = 2Aa{Ri)laiR,)/a{R^)]Aa{Ri), (10) 
w^^Ri) = 0.7Aa{Ri)/a{Ri), (11) 
w^^{Ri) = 0.95aA/i4. (12) 

The coefficients are representative in the range of values 
taken by the observed error bars and the measured hi. In the 
presently used fitting procedure, the additional dependence 
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on h4 and the respective "other" error bar is neglected - the 
tests below show that this leads to satislactory results. Since 
to first order the S2-nioment thus measures the shift from 
a to CT, the parameter a in eq. (13) will be near a ~ 2. We 
have fixed a by requiring that the distributions of Xa- and 
Xh4 have equal width (see below); this results in a = 1.7. 

Finally, we assume that the DF underlying the observed 
kinematics is smooth to the degree that is compatible with 
the measured data. Clearly, unless such an assumption is 
made, it is impossible to determine a function of two vari- 
ables, f{E,L'^), from a small number of data points with 
real error bars. One way to ensure that the DF is smooth is 
to use only a small number of terms in the expansion (5). 
However, this is not a good way of smoothing as it biases 
the recovered DF towards the functional forms of the few 
fk that are used in the sum. A better way of smoothing is 
the method of regularization, as recently discussed by Mer- 
rill (1993) in a similar context. Regularization has tradition 
in other branches of science, and different variants exist; for 
references see Press et al. (1986) and Merritt's paper (1993). 
In the algorithm here, regularization is implemented by tak- 
ing the number of basis functions large enough (typically, 
16 — 20) that the data can be modelled in some detail, and 
then constraining the second derivatives of the composite 
DF to be small. That is we also seek to minimize 



wl{E) X 



(13) 



on a grid of points {Ei,Xj) in the energy and circularity in- 
tegrals. Here we normalize to the isotropic DF to ensure that 
fluctuations in the composite DF are penalized equally at all 
energies; i.e., Wr{E) — l/fiso{E). The constant De is pro- 
portional to the range in potential energy in the respective 
model. 

To find a regularized spherical DF for given kinematic 
data in a specified potential, we thus minimize the quan- 
tity 




+ A^A(/),, 



(14) 



for given regularization parameter A, subject to the equality 
and inequality constraints (6) and (7). For the actual nu- 
merical solution we have used the NETLIB routine LSEI by 
Hanson and Haskell (1981). Once the best model is found, 
we redetermine its quality by evaluating its deviations from 
the actually measured data: 



1 



2 

Xh4 



= 7 



1 



/{AafiR.), (15) 



/(A/i4)'(i?0- (16) 



The parameters a'^-^^Ri) and h\l\Ri) are determined by 
fitting a Gauss-Hermite series to the velocity profiles of the 
best fitting model; I is the number of kinematic data points. 



4.2 Tests with model data 

We have tested this method by applying it to kinematic 
datasets generated in the following Monte Carlo-like way. 
First, velocity dispersions and /14-parameters are calculated 
from a theoretical DF of specified anisotropy in a known 
potential, and are interpolated to the radii Ri/r,j at which 
observed data points are available for NGC 6703. Error bars 
at these Ri are taken to be either the measured error bars 
for NGC 6703 (a realistic dataset), or are taken to be inde- 
pendent of radius with Act = 3 kms~^ and A/i4 — 0.01 (an 
idealized dataset). Then pseudo data are generated from the 
model values of a and hi at Ri by adding Gaussian random 
variates with l-a dispersion corresponding to the respec- 
tive Act or A/14 error bar at this point. Figs. ^, ^ show 
datasets generated in this way from a radially anisotropic 
model in a potential of a self-consistent Jaffe sphere with 
rj = 46. "5, GMjjrj = 272kms"\ and a dark halo with 
Wo = 220kms~^, = 56". This potential was chosen be- 
cause it lies in the middle of the range of acceptable poten- 
tials for NGC 6703 (see Sect. 5.2), so that any systematic 
errors in our analysis will be similar for this model and for 
the galaxy itself. 

We have used the regularized inversion algorithm de- 
scribed in the last subsection to analyse several such pseudo 
data sets, and have determined composite DFs as a func- 
tion of the regularization parameter A. The algorithm was 
given 16 basis function models. These included the isotropic 
model and a variety of tangentially anisotropic models with 
different anisotropy radii and circularity functions but not, 
of course, the radially anisotropic model from which the data 
are drawn (see Fig. H and Appendix A). For each composite 
model returned by the algorithm we have determined two 
diagnostic quantities. The first is the mean x^ per a and /14 
data point, Xa+hA^ which measures the level at which this 
model fits the data from which it was derived. The second 
is the rms deviation between the returned DF and the true 
DF of the model from which the data were drawn, in some 
specified energy range. 

The use of these diagnostic quantities requires some fur- 
ther comments. As usual, the number of degrees of free- 
dom in such a regularized inversion problem is not well- 
determined. For near-zero A, in the case at hand we can 
adjust 16 coefficients at and an overall mass scale. The to- 
tal number of a and ft4-data points is 70; thus in this limit 
the number of degrees of freedom is 53. For large A, on the 
other hand, the recovered DF will be linear in both E and 
X and the values of the at are essentially fixed. Then the 
number of degrees of freedom approaches 69. In both cases, 
it is of the order of the number of data points. Hence the 
use of x^+h4 s-nd /14 data point instead of a reduced 

The second diagnostic measuring the accuracy of the re- 
covered DF must clearly depend on the range in energy over 
which it is calculated. Typically, a kinematic measurement 
at projected radius Ri contains information about the DF 
in a range of energy above the energy of the circular orbit 
at radius Ri. The precise upper end of this range is model- 
specific; it depends on the potential, the kinematic proper- 
ties of the DF itself, and, through the projection process, 
also on the stellar density profile. The use of the outermost 
kinematic data points thus contains, explicitly or implic- 
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Figure 4. Model analysis of pseudo-data generated from a radi- 
ally anisotropic model DF with the sampling and the measured 
error bars of NGC 6703 (see text). The full curves show the true 
profiles of projected velocity dispersion, line-of-sight velocity dis- 
tribution parameter h^, and anisotropy parameter /3 of the un- 
derlying DF. The dashed and dotted lines show the cr, /i4, and /3 
profiles of five regularized composite DFs which were computed 
by the method of Sect. 4.1 with A = 10~*. One of the four curves 
in each panel corresponds to the data points actually shown in 
this figure; the other three curves derive from statistically identi- 
cal kinematic data sets with randomly different values for a and 
h4 within the same (Gaussian) errors. 



itly, assumptions on the radial smoothness of these quan- 
tities. In a typical elliptical galaxy problem, the DF must 
be known fairly accurately at around 3Ri for the projected 
kinematics at Ri to be securely predicted, and for radially 
anisotropic models, even the values of the DF near lOiii can 
make some difference. The RMS residuals in the recovered 
DF given below have therefore been calculated in a range of 
energies extending from "l>(r = 0) to ^{r — 3R,n), where the 
last kinematic data point in the data sets used is located at 
i?™. = 1.68rj. 

Figure ^ shows these quantities as a function of the regu- 
larization parameter A for the two pseudo data sets shown in 
Figs. ^, |5| For small values of A, the composite models fit the 
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Figure 5. Model analysis of pseudo-data generated from the 
same radially anisotropic model DF as in Fig. ji] but assuming 
radially constant error bars Act = 3 km s~^ and A/i4 = 0.01. The 
full curves again show the true ct, /i4, and /3 profiles of the under- 
lying DF. The dashed lines show the projected ct and /14 profiles 
derived by the method of Sect. 4.1 from the pseudo-data shown, 
with A = 6x 10~^. The results obtained for these projected quanti- 
ties from statistically identical data sets are now nearly identical. 
The dotted curves show the uncertainties that remain in the de- 
projected quantity [3 even with such small error bars. 



data accurately, but the recovered DF is not very accurate 
because it contains large spurious oscillations depending on 
the particular values of the data points. For large values of 
A, the models are so heavily smoothed that they neither fit 
the data well, nor do they represent a good approximation 
to the true DF. The optimal regime is where the smooth- 
ing is large enough to damp out the spurious oscillations, 
but still permits resolving the important structures in the 
underlying DF. For values of A in this regime the fit to the 
data is still satisfactory, and the representation of the DF is 
optimal. Fig. ^ shows that the RMS residuals of the DF go 
through minima at values of A ~ 2 x 10~^ for the pseudo 
data in Fig. ^ and A ~ 5 x 10~^ for those in Fig. |5| [when 
De = 1; see below eq. (13)]. 

The shapes of the Xa+hAi^) curves were found to be al- 
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Figure 6. Model results as a function of the regularization pa- 
rameter A, for the two pseudo-data sets shown in Figs. |4|and|5[ 
The two data sets are flagged by triangles (Fig. ^ and pentagons 
(Fig. ^). The upper two curves marked by the open symbols show 
the total P^r o" and h4 data point of regularized composite 
models with 16 basis functions. The lower two curves marked by 
the filled symbols show the rms deviation between the recovered 
composite model DFs and the true DP from which the data sets 
were generated. This rms deviation was evaluated on a grid in 
energy and angular momentum corresponding to radii < 3 times 
the radius of the outermost data point. 



ways similar to the upper curves in Fig. |^. The shapes of the 
corresponding lower curves in Fig. |^ are more variable. The 
resulting optimal values for A can vary, depending on the 
random realisation of the data within the assumed Gaussian 
errors, as well as on the distribution function and potential 
from which the model values are drawn. We have therefore 
investigated 100 realisations of data generated from each of 
several model DFs and potentials. Based on these experi- 
ments we have fixed optimal values of A = 1 x 10~^ and 
A = 6 X 10~® for data with error bars like those in Figs. ^ 
and ^ respectively, for all models that include a dark halo 
component. 

We have done similar experiments with a self-consistent 
model that is radially anisotropic near the center and tan- 
gentially anisotropic in its outer parts, such as might be 
relevant in tests for dark matter at large radii. In these tests 
we have used 20 basis functions. To match the corresponding 
pseudo-data with similar accuracy requires smaller optimal 
values of A (~ 3 X 10~^), so as to compensate for the larger 
derivatives in eq. (13). 

Fig. ^ shows the cumulative distributions of x'i+h4 ^'iid 
of the RMS residual in the DF, as described, for the dynami- 
cal models recovered from several such sets of Monte Carlo 
data. From the top panel it is seen that our fitting algo- 
rithm will match kinematic data in a known potential with 
xl+h4, < 1 about 60-70% of the time, and with xl+h4, > 1-28 
only < 5% of the time. These numbers are similar to those 
expected from a ^^-distribution with 70 degrees of freedom 
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Figure 7. Top: The cumulative distribution of the total normal- 
ized XCT+h4' ^^'^ dynamical models recovered from 100 random 
Gaussian data sets derived from the underlying true DF and the 
observational errors in Figs. ^ (solid line) and ^ (dashed). Bot- 
tom: Cumulative distribution of residuals between true and re- 
covered distribution function, evaluated on a grid extending to 
three times the radius of the last data point, for the same 100 
datasets from both models. Also shown in the top panel are the 
cumulative x^_|_jj4^distributions for data drawn from radially an- 
isotropic models in the two potentials corresponding to two of 
the extreme solid lines in Fig. |ll| with the same error bars as 
in Fig. ^ (dot— dashed lines), and for a self— consistent model with 
more complicated anisotropy structure (dotted line). 



(for Gaussian data), which should describe the statistical 
deviations of the Monte Carlo data points from the under- 
lying true DF. If in modelling the data for NGC 6703 a level 
of XCT-i-h4< 1-28 cannot be reached, the assumed potential is 
not correct with 95% confidence. 

The bottom panel of Fig. ^ shows the distribution of 
residuals in the recovered DF for one model. There are three 
factors which limit the degree to which a given underlying 
DF can be recovered. The first is determined by the data, 
i.e., the size of the error bars, the sampling, the fact that 
there are no measurements at large radii, etc. The second 
is the level of detail that can be resolved by the modelling, 
given the finite number and the particular form of the basis 
functions. The third is the amount of small-scale gradients 
in the model itself. Figs. ^ and |^ show that, if the potential is 
known, data like those for NGC 6703 (Fig. allow recover- 
ing reasonably smooth DFs to an RMS level of < 12% inside 
3i?m about 60 — 70% of the time. Data with much smaller 
error bars but the same sampling (Fig. |^ would give an RMS 
level of < 10% (< 8% for the other two halo models shown 
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in the top panel of Fig. |7|). In the radial-tangential self- 
consistent model the DF to 3-Rm is recovered with an RMS 
accuracy of 16% and 35% for data with error bars like those 
in Figs. ^ and 0, respectively. These comparisons and the 
results of Fig. Mshow that the former values are dominated 
by the measurement uncertainties rather than the resolu- 
tion in the modelling. This conclusion would be different for 
highly corrugated true DFs. However, we would not be able 
to recover such DFs from realistic data in any case. 

The kinematics of the regularized composite DFs de- 
rived from pseudo-data with the chosen optimal A values 
are shown in Figs. ^ and |^. (For brevity, such model DFs 
obtained with near-optimal A will henceforth be denoted as 
"best-estimate models" .) The good match to the data points 
is apparent. The differences in the intrinsic anisotropy pa- 
rameter P between the recovered models and the true model 
are larger. For the data set generated with the observed er- 
ror bars of NGC 6703 the recovered P values are uncertain 
by A/3 ~ ±0.1, and slightly more at the largest radii where 
the observational errors are large. For the pseudo-data with 
the small error bars in Fig. ^ this uncertainty is reduced. 

4.3 Constraining $ 

So far we have shown how the stellar DF can be recovered 
from VP-data in a known spherical potential. In this Section 
we return to the discussion of Section 3 and investigate the 
degree to which the gravitational potential itself can be con- 
strained. Eventually, it will clearly be important to answer 
the theoretical question of whether, in principle, the gravita- 
tional potential is nearly or even uniquely determined from 
the projected DF A''(r, However, we shall not attempt 
to do this here; previous theoretical work suggests that the 
range of potentials consistent with ideal data is small (G93, 
Merritt 1993, see Section 3). More relevant at the moment 
perhaps is the more practical question of how well the po- 
tential can be constrained from real observational data with 
realistic error bars, finite radial extent, and limited sam- 
pling. What is the range of pairs (/, $) that correspond to 
the same data? How does this range shrink as the data im- 
prove? 

Here we investigate these questions with a view to the 
analysis of the NGC 6703 data below. We again use the 
model underlying Figs. ^ this is a radially anisotropic DF 
in the potential of a self-consistent Jaffe sphere and a dark 
halo with parameters GMj/rj — 272kms~^, rc/rj = 1.2 
and Vo/{GMj/rj) = 0.808. As before, random Gaussian 
data sets were generated from this model, with the positions 
and error bars of the data points (i) as in Fig. ^ (ii) as in 
Fig. ^. The last data point is at R,n = 1.68rj, as for the 
NGC 6703 data, well beyond two effective radii. However, 
contrary to Section 4.2, the data sets used in the following 
were specially selected such that Xct — 1 and Xm — 1- (It 
turns out that the data points shown in Fig. ^ have less than 
5% probability according to Fig. ^). 

For these pseudo data we have determined best- 
estimate DFs in a number of assumed potentials, includ- 
ing the underlying true potential. The potentials were cho- 
sen such that (i) they correspond to approximately constant 
mass-to-light ratio for r <C rj, and (ii) they form a sequence 
of varying true circular velocity Vc{Rm) at the radius of the 
last data point. The sequence, with the slightly different ve- 
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Figure 8. Goodness of fit X^^f^4 foi' best— estimate DFs derived 
in a sequence of luminous plus dark matter potentials. The se- 
quence of potentials, with rotation curves similar to those shown 
in Fig. is here parametrized by the total circular velocity at 
the last observational radius, Vc{Rm)- The DFs were fitted to 
pseudo data generated from a model with Vc{Rm) = 242kms~^. 
Solid pentagons: data points with error bars as in Fig. ^. Open 
pentagons: data points with error bars as in Fig. ^. 

locitynormalisations appropriate for NGC 6703, is shown in 
Fig. hll. For each potential, specified by the selected values 
of the parameters rc/rj and v'q/(GMj /rj), we determined 
the goodness-of-fit Xct+/i4 ^ ^ function of the velocity scale 
GMj /rj from the corresponding best-estimate models. Fi- 
nally, we computed the x^^+hi of best-estimate DF for 
the optimum velocity scale. This optimal velocity scale usu- 
ally turns out slightly different for the two pseudo data sets. 

Fig. ^ shows the goodness-of-fit Xa+hi of the opti- 
mal sequence of best-estimate models as a function of 
Vc{Rm), as determined for both model datasets. The un- 
derlying true potential has Vc{Rm) = 242 km s^^ with the 
adopted parameters. Fig. ^ shows that only potentials with 
Vc{Rm) < 230kms~^ can be ruled out (have less than 5% 
probability according to Fig. ^) from the data with "re- 
alistic" error bars. This conclusion is not surprising given 
how large these error bars are. However, the surprise is that 
also with "idealized" data there remains a range of poten- 
tials with larger Vc(Rm) than the true value, in which these 
data can be fit no worse than in the true potential: formally, 
235kms~^ < Vc{Rm) < 285kms~^ with 95% confidence. 

Fig. ^ shows the predicted kinematics and the intrinsic 
anisotropy of the four models in Fig. ^ that match the "ideal- 
ized" data set with Xa+h4 — 1- The best-estimate DFs in the 
potentials with higher Vc{Rm) than that of the true potential 
achieve their good fit to the data by the following means; 
Inside rj they compensate the potential's higher circular 
speed by a larger radial anisotropy, which leads to slightly 
larger /i4 -values (cf. Section 3). This effect is too small to 
be detected even with the small error bars. Outside rj, they 
compensate by a radially increasing tangential anisotropy. In 
this way, the velocity dispersions near the edge of the model 
can be lowered, because the number of high-energy orbits 
coming in from outside is reduced. Such orbits would con- 
tribute large line-of-sight velocities near their pericentres. 
The hi values in the region concerned are also only barely 
affected, because the extra tangential orbits near radius r 
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Figure 9. Model fits to pseudo velocity dispersion and h4 data 
with "idealized" error bars (see text, and Fig. The full line 
shows the kinematics of the model from which these pseudo data 
were generated. The dashed lines show the kinematics of best- 
estimate DFs in a sequence of potentials with increasing halo 
circular velocity, all of which fit the data perfectly with X^+/i4 — 
1. The intrinsic anisotropy of these models (dotted lines) and of 
the true model are shown in the bottom panel. 



and the lack of higher energy radial orbits from beyond r 
nearly compensate. 

Clearly, this mechanism will work less well both when 
the difference between the true and attempted circular ve- 
locity curves increases, and when the radial extent, sam- 
pling, and quality of the data points improve. With data 
yet better and more extended than those in Fig. ^, some 
of the potentials consistent with the presently used model 
data could probably be ruled out. Thus it appears that, if 
the asymptotic circular velocity is constant, then the value of 
that constant can be determined accurately with sufficiently 
high quality data. 

On the other hand, the potentials we have used are 
very simply parametrized functions. There might well exist 
more complicated potential functions, whose circular veloc- 
ity curves differ in only a restricted range of radii, that would 
be impossible to distinguish even with extremely good data. 



To test this, we have constructed a truly idealized data set 
of two times seventy data points with "idealized" small error 
bars as in Fig. ^ evenly spaced in radius, and extending to 
6rj. A model with potential differing only in the halo core 
radius (66% of the true value) was found to fit even these 
data with Xa+hA = 0.97, while for a model with different 
halo core radius and different asymptotic circular velocity 
(by 30kms~^) a satisfactory fit could not be found. 

We draw the following conclusions from these experi- 
ments: 

(i) Velocity profile data with presently achievable er- 
ror bars contain useful information on the gravitational po- 
tentials of elliptical galaxies. In particular, const ant-M/L 
models are relatively easy to rule out once the data extend 
beyond 2Re. The examples that we have studied in detail, 
tuned to the NGC 6703 data, certainly belong to the less 
favourable cases, because the dispersion profile is falling. 

(ii) The detailed form of the true circular velocity curve 
is much harder to determine. Conspiracies in the DF are 
possible that minimize the measurable changes in the line 
profile parameters. A good way to parametrize the results is 
in terms of the circular velocity at the radius of the outer- 
most data point. With presently available data, this can be 
determined to a precision of about ±50kms~^. 

(iii) Better results can be expected from higher-quality 
data, of the sort one could expect from the new class of 10m 
telescopes. However, even with such data some uncertainty 
on the detailed circular velocity curve will remain - regard- 
less of whether or not in theory the potential is uniquely 
determined from the projected DF N{r,Vii). Therefore, the 
combination of the type of analysis presented here with other 
information (e.g., from X-ray data) will give the most pow- 
erful results. 



5 THE ANISOTROPY AND MASS 
DISTRIBUTION OF NGC 6703 

5.1 Constant— M/L model fits 

The SB-profile of NGC 6703 is well fitted by a Jaffe mo- 
del. The largest local residuals are < 15% around R ~ 25", 
< 10% around ii ~ 40", and smaller elsewhere, in particular 
at i? > 60" . A non-parametric inversion of the the surface 
brightness profile showed that the deviations from a Jaffe 
density law are not significantly larger than those quoted 
in SB. In the curve of growth, measuring the total luminos- 
ity inside R, the residuals are everywhere less than ~ 2% 
(Fig. 1^). Since the potential "l>(r) is determined by the to- 
tal mass Af(< r), we can thus to good approximation use a 
Jaffe model for the gravitational potential of the stars. This 
enables us to also compare our kinematic data with a large 
set of self-consistent dynamical models from Jeske (1995). 

We have fitted all models from this database to the 
observational data for NGC 6703, taking the fitted Jaffe 
radius rj = 46. "5. All data points from Fig. ^ were included 
and weighted equally with their individual error bars, except 
for two adjustments, (i) The error bars of the three /i4-points 
near 25" have been set equal to their standard deviation, and 
(ii) the error bars of the innermost eleven /14-points have 
been set to twice the measured values. These modifications 
prevent the total xL to be dominated by these points which 
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Figure 10. x|^-XCT"diagram from fitting dynamical models to 
tfie kinematic data for NGC 6703. Tlie crosses show fit results for 
a variety of anisotropic, self-consistent models of a Jaffa sphere, 
whose density profile is a good approximation for the luminous 
matter in NGC 6703. They fall upwards and to the right of a 
curved envelope that separates them clearly from a perfect fit, 
showing that no self-consistent model can simultaneously fit both 
the dispersion and line profile shape data. The heavy dot is the 
best-estimate self-consistent model in the stars— only potential, 
obtained with the method of Section 4. The squares show a num- 
ber of dynamical models with a dark halo; these are consistent 
with the data. 

are of no consequence for the halo of NGC 6703. Moreover, 
systematic effects may play a role for the innermost data 
points (Sect. 2). 

The velocity scale of each model is matched to opti- 
mize the fit to the cr-profile; that to the /14-profile is already 
fixed with no extra assumption. The resulting values of , 
Xct are normalized by the number of fitted data points. In 
Fig. ^ we plot a Xh^-XCT-diagram for all self-consistent mo- 
dels from Jeske (1995). The values are again normalized 
to the number of data points. Also plotted on the figure are 
best-estimate models in a number of potentials, constructed 
with the technique described in Section 4. The optimized 
self-consistent model (the heavy dot) lies on the bound- 
ing envelope of the self-consistent models; it has Xct = 2.6, 

= 1.3. The squares show the normalized x^-values for 
several models with a dark halo whose rotation curves are 
shown by the full lines in Fig. 9. For most of these Xct — 1-3, 

- 0-8. 

Fig. hd shows that no self-consistent model will Rt the 
data: all self-consistent models are clearly separated by a 
curved envelope from the lower left hand corner of the dia- 
gram, which corresponds to a perfect fit. Either the velocity 
dispersion profile is matched reasonably well, but then the 
line profiles cannot be reproduced, or, when the /14-profile 
is fitted accurately, the dispersion profile is poorly matched. 
The cure for the discrepancy is to raise both a and /14 at 
large radii. Thus, according to Sect. 3 above, we require ex- 
tra mass at large r. NGC 6703 must have a dark halo. 

5.2 Dynamical models with dark halo 

We will now derive constraints on the gravitational poten- 
tial of NGC 6703 within the framework of the parametric 
mass model of eqs. (2)-(4). In doing this we have in mind 




Figure 11. Rotation curves for a sequence of gravitational po- 
tentials (stars plus dark halo) used in the analysis of NGC 6703. 
The full lines show rotation curves that are consistent with the 
NGC 6703 kinematic data inside the 95% confidence boundary at 
A = 10~^ (open symbols in Fig. hd). The other line styles show 
rotation curves inconsistent with the data; among these is the 
constant-M/L model with no dark halo (short-dashed). 



the following working picture: The stellar component is as- 
signed a constant mass-to- light ratio T, chosen maximally 
such that the stars contribute as much mass in the center 
as is consistent with the kinematic data. The model for the 
halo incorporates a constant density core, and its parame- 
ters are chosen such the halo adds mass mainly in the outer 
parts of the galaxy if that is necessary. This is similar to 
the maximum disk hypothesis in the analysis of disk galaxy 
rotation curves. Within this framework we can determine 
the maximum stellar mass-to-light ratio, ask whether it is 
reasonable, and constrain the halo parameters. 

Because determining the two potential parameters Vc 
and Wo together with the model velocity scale (equivalent 
to the stellar mass-to-light ratio) is a three-dimensional 
problem, we will proceed in steps. Fig. ^ shows circular 
rotation curves for the first sequence of mass models for 
NGC 6703 that we have investigated. In all these models the 
halo contribution becomes significant only outside 20" < Re- 
The corresponding halo core radii mostly lie between 1.2rj 
and 1.7rj, i.e., 55" and 80". These values are relatively large 
because of the falling dispersion curve in NGC 6703. This 
implies that we can reliably determine only one of the halo 
parameters for this galaxy. 

We have chosen the circular velocity Vc{Rm.) at the ra- 
dius Rm of the last kinematic data point as this parameter. 
The sequence in Fig. |ll] was constructed such as to vary 
Vc{Rm) and the rotation curve outside ~ Re, while leaving 
the central rotation curve nearly invariant. However, when 
the model velocity scales are optimized in the determina- 
tion of the best-estimate DFs, the optimal velocity scale 
is found to correlate inversely with Vc{Rm)- The rotation 
curves in Fig. |l^ are plotted with their optimal scaling; then 
the central rotation curves are no longer identical, but be- 
come scaled versions of each other. 

In each of the corresponding potentials, we have con- 
structed the best-estimate DF with optimal velocity scaling, 
as described in Section 4. We first fit a composite DF to the 
velocity dispersion and line profile shape data for a series of 
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Figure 12. Quality with which the kinematics of NGC 6703 
can be fitted in different potentials. The figure shows the av- 
erage psr o"- and data point, of the best— estimate distri- 
bution function fitted to the velocity dispersion and line profile 
data, as a function of the assumed potential's circular rotation 
velocity at the observed radius of the last kinematic data point. 
Filled symbols show best-estimate models derived with the opti- 
mal A = 10~* determined in Section 4.2; open symbols represent 
models derived with A = 10~^. The self— consistent (A//L = const.) 
and the Vc = const, models are marked separately. The horizon- 
tal dashed line shows the 95% confidence boundary derived from 
Fig. |. 



values of the unknown velocity scale. From this sequence of 
models, we determine the optimal value of the model's ve- 
locity scale in this potential, and then recompute the best- 
estimate model with this velocity scale. In the following, 
when speaking of the best-estimate DF for a given poten- 
tial, we will always imply that the velocity scale has been 
optimized in this way. In the fitting procedure we have used 
a regularization parameter A = 10"^; this was found appro- 
priate in Section 4.2 for the error bars and sampling of the 
NGC 6703 data. Compared with the tests in Section 4, we 
have included a few extra basis functions (total K — 20) to 
resolve the (possibly not real, cf. Section 2) high-frequency 
structure in the center of NGC 6703. 

Fig. ^ shows, as a function of Vc{Rm), the average 
per a- and hi data point of the respective DFs so obtained. 
The connected solid symbols remesent the sequence of po- 
tentials corresponding to Fig. hll. The potential with con- 
stant mass-to-light ratio appears in the upper left-hand 
corner in Fig. [12| it is inconsistent with the data by a large 
margin even for optimum velocity scale (see Fig. |^. The 
best-fitting potential with a completely flat rotation curve 
has an optimal value of Uc = Vc{Rm) = 254kms~^ and 
Xa+hi — 1-17. Thus it does not provide the best possible fit 
but is consistent with the data. Of the stars plus dark halo 
models illustrated in Fig. |ll|, those models in the sequence 
with Vc(Rm) = 210 - 285 kms"^ have xl+hi < 1-28. This 
is consistent with the results of Section 4.3, from which we 
would expect that the NGC 6703 data can be fit by a range 
of gravitational potentials. Models in the sequence outside 
this range of Vc{Rm) are inconsistent with the kinematic 
data at the 95% confidence level (cf. Fig. 

In a second step we have analyzed a more complete set 
of potentials in a suitable part of the (uo, rc)-plane (Fig. p^ . 




Figure 13. The {vo,rc) halo parameter plane. Values of Vc are 
scaled with respect to rj = 46". 5. The luminous plus dark matter 
models investigated are shown as plus signs. The contours show 
lines of constant X^+hi obtained by interpolating between the 
model values. Acceptable potentials lie in a band extending from 
low vo and low rc to high vq and high re- Models at the upper 
left are ruled out because they do not contain enough mass at 
large radii. Models at the lower right are ruled out because no 
satisfactory fit can be found for any constant value of the stellar 
mass— to-light ratio. 



The fitting results in these potentials are shown by the iso- 
lated filled symbols in Fig. |l^. Combining with the previous 
results these allow us to investigate the range of accept- 
able potentials for NGC 6703 more fully than with the one- 
dimensional sequence in Fig. |ll| 

Fig. shows contours of constant Xa+h4 in the {vo,rc)- 
plane. The most probable potentials lie in a band extending 
from low Vo and low rc to high no and high rc. Thus, as al- 
ready discussed above, it is not possible to determine both 
halo parameters in the NGC 6703 case. However, potentials 
in the band of most probable («o, ^c) all have circular veloc- 
ities Vc{Rm) in the same range of 250 ±40 km s~^ as before. 

In fact, the best-fitting velocity scales of all models in 
Fig. ^ turn out such that the resulting values of Vc{Rm,) are 
in the range [189, 318] km s~^. The fitting procedure tends 
to move Vc{Rm) into the correct range even when no satis- 
factory DF can be found. E.g., some models in the lower 
right of Fig. ^ with x1+h4 — 1-5 appear in Fig. |l^ at 
Vc{R,n) = 250 - 300kms"\ some with xI+m Z 2.4 (not 
shown) at Vc{Rm) = 280 — 340kms~^. All models shown in 
the upper left of Fig. |l^ fall near the line defined by the se- 
quence of models discussed at the beginning of this section. 
The fact that Vc{Rm) varies relatively little for a wide range 
of luminous matter plus dark halo models suggests that our 
results, when expressed in terms of this parameter, are not 
sensitive to the choice of halo model in eqs. (3), (4). 

From Fig. |l^ we conclude that the true circular ve- 
locity of NGC 6703 at 78" is Wc(i?m) = 250 ± 40kms-^ 
(the formal 95% confidence interval obtained from the filled 
symbols, according to Fig. |^. In Section 4.3 we found, 
however, that most of this indeterminacy is towards large 
circular velocities, at least for a galaxy with a dispersion 
curve like that of NGC 6703. By contrast, in potentials with 
lower values of Vc{Rm) than the true circular velocity, it 
quickly becomes impossible to find a satisfactory DF. Based 
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on these results, the lower values in the quoted range of 
Vc{Rm) ~ 250 ± 40kms~^ appear to be the more probable 
ones. 

The open symbols in Fig. ^ show that the range of 
potentials consistent with the data is enlarged only slightly 
when the DF is allowed to be less smooth. These points are 
from best-estimate models derived in the sequence of po- 
tentials of Fig. |ll| with A = 10~^ instead of the optimal 
A = 10"". The resuhing curve which connects the A = 10 
models in Fig. |l^ surrounds the corresponding A — 10~" 
curve. Generally it appeared that, in the models obtained 
with A — 10~^, the DF came close to zero more easily and 
more often. The last two facts, when taken together, sug- 
gest that, in addition to the data themselves, the positivity 
constraints on the DF play an important role in determining 
the boundary of the region in (/, $) space that is consistent 
with given kinematic data. 

Fig. |l^ shows best-estimate models in some of the 
potentials consistent with the kinematic data (at A = 
10""). The three full lines are models with Vc{R,n) = 
(218, 231, 242) kms"\ the two dotted lines are models with 
Vc{Rm) = (253, 277) kms"'^. Also shown is the best- 
estimate model in the potential of only the stars with con- 
stant M/L (long-dashed), and that in the best potential 
with Vc — const. = 254 km s"'^ (short-dashed). 

The dispersion profile of the best-estimate self- 
consistent model falls clearly below the data at both inter- 
mediate and large radii, as does its /14-profile. From the dis- 
cussion in Section 3, this is a clear sign of extra mass at large 
radii. The models including dark halo contributions mainly 
differ in the outermost parts of the velocity dispersion pro- 
file. As expected, those with the highest velocity dispersions 
at large radii correspond to the potentials with the largest 
asymptotic circular velocities. This again suggests that with 
smaller error bars at large radii and with spatially more ex- 
tended data, we will be able to significantly narrow down the 
uncertainties in the halo parameters. The model with every- 
where constant rotation speed is constrained tightly by the 
kinematic data in the central parts, where the p (x~ pro- 
file is presumably dominated by the stars. It then has some 
difficulties both with the /14-values at intermediate radii and 
the velocity dispersions at large radii. 

The lower part of Fig. ^ shows that, in order to match 
the observed kinematics in a potential with large circular ve- 
locity, the DF must become rapidly tangentially anisotropic 
at the radii of the last data points and beyond. As discussed 
in Section 4, this is different from the better known efi'ect of 
increasing the projected dispersions in a potential with not 
enough mass at large radii, by making the DF more tangen- 
tially anisotropic. Here, without the tangential anisotropy, 
our models would predict too high values of the velocity 
dispersion, because too much mass at large radii is implied. 
The tangential anisotropy at radii outside those reached by 
the observations reduces the number of orbits that come into 
the observed range, orbits which would contribute relatively 
large line-of-sight velocities near their turning points. In this 
way the velocity dispersions can be reduced down to the ob- 
served values. Clearly, more spatially extended data would 
reduce this freedom. 

From Fig. ^ we conclude that the stellar distribution 
function in NGC 6703 is near-isotropic at the centre and 
then changes to slightly radially anisotropic at intermedi- 
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Figure 14. Dynamical models for the kinematics of NGC 6703 
in several luminous plus dark matter potentials, compared to 
projected velocity dispersion (top panel) and VP-shape param- 
eter /i4 (middle panel). The bottom panel shows the models' in- 
trinsic anisotropy parameter /3(r), with the same linestyles: self- 
consistent model (stars only; long— dashed), Vc = const, model 
(short— dashed) , three models with vdJS") in the lower part of 
the acceptable range (full), and two models with ■Uc(78") in the 
upper part of this range (dotted lines). 



ate radii (/? = 0.3 - 0.4 at 30", /3 = 0.2 - 0.4 at 60"). It is 
not well-constrained near the outer edge of the data, where 
formally {3 = —0.5 1-0.4, depending on the correct poten- 
tial in the allowed range. However, the models with large 
asymptotic halo circular velocities shown in Fig. |l^ appear 
less plausible, because they are the models with the most 
rapidly increasing velocity dispersions outside R ~ 60" . The 
same models also show the most rapid increase in tangen- 
tial anisotropy at and beyond Rm ~ 78", which again ap- 
pears a priori implausible, because it implies rapid changes 
in the DF just outside the observed range. The combined 
signature of both eff'ects is strongly reminiscent of Fig. ^, 
where it was clearly an artifact of the limited radial range 
of the data. If this assumption is correct, one would again 
conclude that loweT-Vc{Rm) models in the formal range are 
favoured. Fig. ^ shows the recovered DF for the potential 
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Figure 15. Distribution function in the energy— circularity plane 
derived for NGC 6703, in the luminous plus dark matter potential 
with Vc{Rm) = 231kms~^ shown as the middle solid line in the 
upper panel of Fig. Energy is specified by the value of r(E), 
where E = 3>(r), in units of the fitted Jaffe radius rj, and the 
DF is given in arbitrary logarithmic units. The last measured 
kinematic data point is located at \gr{E) /r j = 0.23. 



with Vc{Rm) = 231kms~^ in Fig. Note, however, that 
all models shown in Fig. ^ except the self-consistent one 
are formally consistent with the presently available data. 
From the constraints on the circular velocity Vc(Rn 



210- 
is 



-290 kms ^ the range in mass inside Rm = 13.5 h^^ kpc 



M(< Rm) = 1.6-2.6 X 10^ 



(17) 



where /iso = iJo/50km/s/Mpc. The total mass in stars in- 
side this radius is 8 — 9 x 10^° Mq, assuming constant mass- 
to-light ratio T and a maximum stellar mass model, and 
taking an average value from the models consistent with the 
kinematic data. The radial run of the luminous, dark, and 
total mass is shown in Fig. |l6| for the models that span the 
allowed range according to Fig. |l^. After dividing by the lu- 
minosity Lb {r) for the stars, the mass-to-light ratios shown 
in Fig. O result. Between the centre and the last data point 
r = 78"^2.6i?e, the mass-to-light ratio of NGC 6703 rises 
by a factor of 1.6 — 3. 
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Figure 16. Luminous, dark, and total mass as a function of ra- 
dius for the range of acceptable models of NGC 6703, according 
to Fig. |l2| (short dashed, dotted, and dash— dashed or full lines, re- 
spectively). Mass distributions in which a DF with X^^h4 — ^■^'^ 
(including 87% of the distribution in Fig. ^, 1.5<t) are coded by full 
lines, those with X^^hi — ^'"^^ (including 95%, 2. Oct) by dash- 
dashed lines. The vertical line denotes the position of the last 
kinematic data point. At this radius, the luminous mass fails by 
at least a factor of 1.6. 




Figure 17. The B-band mass-to-light ratio of NGC 6703. The 
solid and dash-dashed lines (coding as in Fig. |l^) are derived from 
the dynamical models that span the range of acceptable Vc{Rm) 
in Fig. The central mass— to-light ratio is T = 3.3, that at the 
position of the last kinematic data point at 78" (vertical bar) is 
in the range of T = 5.3 — 10. 



5.3 Uncertainties 

There are a number of possible sources of systematic error 
which would affect the mass-to-light ratio derived for NGC 
6703. Most of the errors so introduced are probably small 
compared to the considerable uncertainty arising from the 
kinematic measurement errors and limited radial sampling, 
discussed above. One systematic error on the absolute mass- 
to-light ratios in NGC 6703 comes from the uncertainty in 
the distance, although this does not change the ratio of outer 
to central values. A further systematic effect on this ratio 
can be introduced by the sky brightness level: If this is in- 
creased by 2 — 3%, the fitted r j decreases to 35" . To the ex- 
tent that the outermost kinematic data point at 78" (which 
then moves to larger r/rj) is in the flat part of the circu- 
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lar rotation curve, the inferred M/L changes only to second 
order because the luminosity inside 78" remains essentially 
unchanged. The same is true if the sky value is decreased by 
1%, in which case the fitted rj increases to 54". 

In the previous analysis, we have ignored a possible 
small rotation in the outer parts of NGC 6703 (perhaps 
~ 20 — 30kms~^ at i? > 50", but the errors are large; cf. 
Fig. |l]). The simplest possible estimate of the effect of this 
rotation on the derived masses is to replace a ~ 105kms~^ 
in this radial range by (cr^ + ~ IfOkms"^. This gives 

a factor of 1.1, neglecting changes in the model structure 
that would result because the central kinematics remain un- 
changed. 

Next we consider the possibility that NGC 6703 may 
contain a face-on extended disk (see de Vaucouleurs, de Vau- 
couleurs and Corwin 1976). From an _R'^^''-law plus disk de- 
composition, we estimate that the contribution of such a 
disk in the region where we model the kinematics could be 
up to 10-20%. In this case we expect the velocity dispersion 
to be decreased and the hi coefficient to become more pos- 
itive where the disk contributes significant light (Dehnen & 
Gerhard 1994; see NGC 4660 as an example in BSC), most 
likely in the outer parts. 

Similarly, it is conceivable that NGC 6703 is in reality 
slightly triaxial and is seen from a special direction so as to 
appear EO. The likelihood of this is the smaller, the more 
triaxial the intrinsic axial ratios are; thus slightly triaxial 
shapes are the most plausible ones. Again this will imply 
some extra loop orbits seen nearly face-on, similarly increa- 
sing /i4 and decreasing a. 

In both cases, we therefore expect the spherical com- 
ponent in NGC 6703 to have lower and larger a than 
the measured values. A similar analysis of such kinematics 
would, according to the discussion in Section 3, lead to a mo- 
del with greater tangential anisotropy at large radii, with 
the mass distribution less affected. Recall that decreasing 
the mass at large r in a spherical model lowers both a and 
hi. 



6 DISCUSSION AND CONCLUSIONS 

This study is part of an observational and theoretical pro- 
gram aimed at understanding the mass distribution and or- 
bital structure in elliptical galaxies. In the following we first 
discuss general results on potential and anisotropy determi- 
nation, and then proceed to the specific case of NGC 6703. 

6.1 Velocity profiles and anisotropy and mass 

The analysis of the VPs of simple dynamical models in Sec- 
tion 3 has broadly confirmed the conclusions of G93. At 
large radii, where the luminosity profile falls rapidly, the VPs 
are dominated by the stars at tangent point. Then radially 
(tangentially) anisotropic DFs can be recognized by more 
peaked (more flat-topped) VPs with more positive (nega- 
tive) hi than for the isotropic case. Increasing (5 at constant 
potential thus lowers o and increases hi. On the other hand, 
increasing the mass of the system at large r at constant 
anisotropy, increases both the projected dispersion and hi. 
This suggests that by modelling o and hi both mass M{r) 
and anisotropy /3(r) can be constrained. 



In practical applications, such an analysis is compli- 
cated by a number of factors. Radial orbits at large radii 
may lead to increased central velocity dispersions and flat- 
topped central VPs (already pointed out by Dejonghe 1987). 
The former effect can be compensated by a decrease in the 
stellar mass-to-light ratio. The latter is independent of this, 
but can be compensated by changes to the distribution func- 
tion in the inner parts of the galaxy (as in a number of cases 
studied in Sections 4, 5). A more serious uncertainty is in- 
troduced by the possibility of significant gradients in the 
orbit population across the radii of interest. For example, a 
population of high energy radial orbits with pericentres in a 
limited radial range may mimic tangential anisotropy there. 
In many cases it will be possible to exclude such a population 
of orbits by its effects on the VPs at exterior radii, i.e., by si- 
multaneously analysing a number of observed VPs. Yet this 
is least possible precisely at the largest observed radii, where 
mass determination is most interesting. Thus this chain of 
argument suggests (correctly, see below) that the largest un- 
certainty in determining masses and anisotropies in ellipti- 
cals from VP-data is the finite radial extent of these data. 

To analyze realistic data we have constructed an algo- 
rithm by which the distribution function and potential of 
a spherical galaxy can be constrained directly from its ob- 
served a and /14-profiles. To assess the significance of the 
results obtained, we have tested the algorithm on Monte 
Carlo-generated data sets tuned to the spatial extent, sam- 
pling, and observational errors as measured for NGC 6703. 
From such data, the present version of the algorithm recov- 
ers a smooth spherical DP, ~ 70% of the time, to an RMS 
level of better than ~ 12% inside three times the radius of 
the outermost kinematic data point. 

We have used this algorithm to study quantitatively 
the degree to which the gravitational potential can be de- 
termined from such data. Our main conclusion is that veloc- 
ity profile data with presently achievable error bars already 
constrain the gravitational potentials of elliptical galaxies 
significantly. In particular, constant-M/L models are rela- 
tively easy to rule out once the data extend beyond 2Re. 
The examples that we have studied in detail, tuned to the 
NGC 6703 data, certainly belong to the less favourable cases, 
because in this galaxy the dispersion profile is falling. 

A good way to parametrize the results is in terms of the 
true circular velocity Vc{Rm) at the radius of the outermost 
data point, Rm. With presently available data, Vc{Rm) can 
be determined to a precision of about ± < 50kms"^ This 
will improve when high-quality data at several Re become 
available, of the kind expected from the new class of 10m 
telescopes. Apart from the fact that smaller error bars will 
decrease the formally allowed range in Vc(Rm), tests show 
that this range often includes high-Wc(-Rm) models which be- 
come rapidly tangentially anisotropic just outside the data 
boundary. These (not very plausible) models can be elimi- 
nated with data extending to larger radii. 

On the other hand, the detailed form of the true circu- 
lar velocity curve is much harder to determine than Vc{Rm). 
Conspiracies in the DF are possible that minimize the mea- 
surable changes in the line profile parameters. Our tests 
showed that two potentials differing by just the value of the 
halo core radius could not be distinguished even with very 
good data out to QRe . Thus some uncertainty will remain in 
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practise, regardless of whether or not in theory the potential 
is uniquely determined from the projected DF N{r,vn). 

A similar picture holds for the related determination 
of the anisotropy of the DF. For the present error bars in 
the data, /3(r) is relatively well-determined out to about 
half the limiting radius of the observations. Near the edge of 
the data, uncertainties can be large depending on the grav- 
itational potential (recall that in a fixed spherical potential, 
the DF is uniquely determined by the complete projected DF 
N{r,V\\)). Again the unknown nature of the orbits beyond 
the last data point has a large part in this uncertainty. 

Because the largest uncertainties in determining masses 
and anisotropics from VPs occur near the outer radial limit 
of these data, the combination of the type of analysis pre- 
sented here with other information (e.g., from X-ray data) 
will be particularly powerful. 

6.2 The dark halo of NGC 6703 

Fig. |l^ shows that no self-consistent model will Rt the 
kinematic data for NGC 6703. Our non-parametric best- 
estimate self-consistent model is inconsistent with the data 
at the > 99%-level (Figs. |l2[ |). With self-consistent mo- 
dels, either the velocity dispersion profile is matched reason- 
ably well, but then the line profiles cannot be reproduced, 
or, when the /i4-profile is fitted accurately, the dispersion 
profile is poorly matched. The cure for the discrepancy is to 
raise both a and at large radii. Thus, as discussed above, 
we require extra mass at large r. NGC 6703 must have a 
dark halo. 

We have next derived constraints on the parameters of 
this halo as follows. The luminous component is assigned a 
constant mass-to-light ratio T, chosen maximally such that 
the stars contribute as much mass in the center as is consis- 
tent with the kinematic data. Our parametric model for the 
halo incorporates a constant density core, and its param- 
eters (core radius and asymptotically constant circular 
velocity vo) are chosen such the halo adds mass mainly in 
the outer parts of the galaxy if that is necessary. We call 
these models maximum stellar mass models (analogous to 
the maximum disk hypothesis in the analysis of disk galaxy 
rotation curves). 

We find that maximum stellar mass models in which 
the luminous component provides nearly all the mass in the 
centre fit the data well. In these models, the total luminous 
mass inside the limiting observational radius Rm = 78" = 
13.5 h^Q kpc is 9 X W^'^'h^g M©, corresponding to a central 
B-band mass-to-light ratio T = 3.3/i5o Mq/Lq. According 
to Worthey's (1994) models, this is a rather low value for the 
stellar population of an elliptical galaxy and would point to a 
relatively low age (5 Gyrs) and/or low metallicity (less than 
solar). However, the galaxy has a color {B — V)o = 0.93 and 
a central line index Mg2 — 0.280 (Faber et al. 1989) which 
are typical for ellipticals of similar velocity dispersion. 

A larger value of Ho could increase the M/L value and 
alleviate the demands on the stellar populations. However, 
the distance used here (36 Mpc) includes a correction for 
the large inferred peculiar velocity of the galaxy. If we had 
used a distance based on the larger radial velocity in the 
CMB frame, our derived M/L would be even lower. It is 
also implausible that the low central value of M/L stems 
from the contribution of a young stellar population in a disk 



component, which we estimate cannot be larger than 20% of 
the total light (see above). Thus we conclude that the dark 
halo in NGC 6703 is unlikely to have higher central densi- 
ties than inferred from our maximum stellar mass models, 
because otherwise the M/L of the stellar component would 
be reduced to implausibly small values. 

In a recent preprint, Rix et al. (1997) have analyzed 
the velocity profiles of the EO galaxy NGC 2434 with a 
linear orbit superposition method. This galaxy provides an 
interesting contrast to NGC 6703 because it has an essen- 
tially flat dispersion profile. Its kinematics are likewise in- 
consistent with a constant-M/L potential, but are well-fit 
by a model with Vc = const. This can be interpreted as a 
maximum stellar mass model in the sense defined above, in 
which the luminous component with maximal T contributes 
most of the mass inside Re. The kinematics of NGC 2434 
are also well-fit by a range of specific, cosmologically moti- 
vated mass models which, if applicable, would imply lower 
T and significant dark mass inside Re. In NGC 6703, a mo- 
del with Ve = const is formally consistent with the present 
data (within 2a), but it is not a very plausible fit at large R 
and requires large anisotropy gradients between 40" and 70" . 
It will be interesting to see whether future studies confirm 
differences between the shapes of the true circular velocity 
curves of elliptical galaxies. 

Because of the falling dispersion curve in NGC 6703, 
we can determine only one of the halo parameters (The halo 
parameters of the most probable potentials lie in a band 
extending from low vq and low Ve to high vq and high rc.) 
However, the circular velocity Vc{Rm) at the data boundary 
is relatively well-determined for all these models. Thus we 
find (Fig. hd) that the true circular velocity of NGC 6703 at 
78" is Ve{Rn^) = 250 ± 40 km s"^ (formal 95% confidence in- 
terval). Tests on pseudo data have shown that this range of- 
ten includes high-Wc (i?m) models which become rapidly tan- 
gentially anisotropic just outside the data boundary. Such 
models may not be very plausible, so the lower values in the 
quoted range of Vc{Rm) = 250±40kms~^ may be the more 
probable ones. 

Thus, at Rm = 78" = 13.5 h'^Q kpc the total mass en- 
closed is M{< Rm) = 1.6 - 2.6 X W^^h'g Mq, and the inte- 
grated mass-to-light ratio out to this radius is T = 5.3— 10, 
corresponding to a rise from the center by at least a factor 
of 1.6. We have already noted that NGC 6703 is an un- 
favourable case because of its falling dispersion curve. The 
fact that relatively small variations in T can nonetheless 
be detected shows the power of the method. Note that a 
scheme based on the analysis of the line of sight velocity 
dispersions alone (Binney, Davies, Illingworth 1990, van der 
Marel 1991) would conclude that constant mass-to-light ra- 
tio models can provide good fits. 

The stellar distribution function in NGC 6703 is near- 
isotropic at the centre and then changes to slightly radi- 
ally anisotropic at intermediate radii (/3 = 0.3 — 0.4 at 30", 
= 0.2 — 0.4 at 60"). It is not well-constrained near the outer 
edge of the data, where formally p——0.5 — 1-0.4, depending 
on the correct potential in the allowed range. Models near 
the lower end of this range may be consistent with the data 
only because of the limited radial extent of the measure- 
ments. 
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6.3 Conclusions 

In summary, we have shown that the mass distribution M(r) 
and anisotropy structure /9(r) for spherical galaxies can both 
be constrained from VP and velocity dispersion measure- 
ments. NGC 6703 must have a dark halo, contributing about 
equal mass at 2.6Re as do the stars. The circular velocity at 
the last kinematic data point (78" ) is constrained to lie in the 
range 250 ±40 km s~^ at 95% confidence. The anisotropy of 
the stellar orbits changes from near-isotropic at the center to 
slightly radially anisotropic at intermediate radii, and may 
be either radially or tangentially anisotropic at 78". With 
more extended and more accurate data it will be possible to 
considerably narrow down these uncertainties. 

If the results for this gala^xy are typical, they suggest 
that also in elliptical galaxies the stellar mass dominates at 
small radii, and the dark matter begins to dominate at radii 
around 10 kpc. It is important to obtain extended kinematic 
data and do a similar analysis for a number of elliptical ga- 
laxies. When we know the systematics and the spread in the 
circular velocity curves and anisotropy profiles for a sam- 
ple of ellipticals, we will have an important new means for 
testing the currently popular formation theories. 
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APPENDIX A: LIBRARY OF ANISOTROPIC 
SPHERICAL DISTRIBUTION FUNCTIONS 

To understand the connection between anisotropy struc- 
ture and observable line profile shapes we have constructed 
a number of spherical distribution functions of the quasi- 
separable form (Gerhard 1991, G91) 



f{E,L)=g{E) h{x) 



(18) 



where the variable x depends on both energy and angular 
momentum: 

L 



(19) 



Lo is an angular momentum constant, or equivalently, an 
anisotropy radius times a characteristic velocity; Lc{E) is 
the angular momentum of the circular orbit at energy E. DFs 
of the form (18) have the following properties: (i) the circu- 
larity function h{x) has the effect of shifting stars between 
orbits of different angular momenta on surfaces of constant 
energy, while g{E) controls the distribution of stars between 
energy surfaces, (ii) For the most bound stars L < Lc{E) <^ 
Lo; thus the model becomes isotropic in the centre unless 
Lo = 0. (iii) For loosely bound stars L ~ Lc{E) ^ Lo, 
i.e., the angular momentum distribution becomes a func- 
tion of circularity L/Lc{E) which is one-to-one related to 
orbital eccentricity. Outside the anisotropy radius, the DF 
(18) therefore corresponds to an energy- independent orbit 
distribution with constant anisotropy, radial or tangential. 

In these models h{x) is an assigned function built into 
the model to achieve the desired anisotropy (orbit distribu- 
tion). Radially biased distribution functions correspond to 
circularity functions h{x) decreasing with x; for example 



h{x) = ha{x) = (l — x^Y 



(20) 
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Figure 18. Anisotropy parameter projected velocity disper- 

sion (jp, and line profile shape parameter /i4 for several families 
of anisotropic DFs. Left: scale-free radially anisotropic (dotted) 
and tangentially anisotropic (dashed lines) in the potential of a 
self-consistent Jaffe sphere. For these models, circularity func- 
tions of the type (20) and (21) were used with different a and 
c, respectively. The isotropic model is shown for reference (solid 
line). Right: Families of models with anisotropy changing from 
radial to tangential (dotted) or from tangential to radial (dashed 
lines). These were constructed from the same circularity functions 
and a weight factor as in eq. (23). 



The family can also be used to construct tangentially aniso- 
tropic models, such as 



he 



c+{l 



[1 



(1 



(21) 



In these tangentially anisotropic models one cannot choose 
h{0) — unless Lo — 0, otherwise the density at r = would 
be zero. Of course, other forms for h{x) axe possible, such 
as Gaussians. 

Given the assigned function h{x), the integral equation 
for p(r) in terms of f{E, L ) is solved for the derived func- 
tion g{E); see G91 and Jeske (1995). Fig. ^ shows line profile 
shape parameters for representative DFs constructed in this 
way. Fig. ^ shows the anisotropy profiles of two sets of tan- 
gentially anisotropic models. Here the density of stars has 
been taken to be that of a Jaffe sphere, and the potential 
in which the stars orbit is either the self-consistent poten- 
tial or one of the mixed stars plus halo potentials used in 
Sections 4, 5. Sequences like that in Fig. M are used as basis 
functions in the non-parametric analysis in Section 5. 

Models whose anisotropy changes from radial to tan- 
gential or vice versa were constructed by linearly combining 
the above circularity functions with energy-dependent co- 
efficients. In this way one obtains DFs of the more general 
form 



f{E,L)=g{E) h{E,x) 



(22) 



For self-consistent Jaffe models we used energy-dependent 
coefficients fJ,{E) of the following form: 



H(E) ~ — -\ arctan 

2 TV 



kE'^ 
2 



E 



£2 



(23) 



The parameters E and k determine the orbital energy near 
which the anisotropy transition occurs, and the width of the 
transition. A similar function fJ,{E) was used for models in 
halo potentials. Figure ^ shows the intrinsic and projected 
properties of a number of DFs of this kind, constructed in 
the self-consistent potential of a Jaffe sphere. Notice the 
wide variety of kinematical profiles that can be constructed 
in this way. 



APPENDIX B: TRANSFORMING TO LINEAR 
KINEMATIC DATA 

In velocity line profile measurements, the depth of an ab- 
sorption line in a spectral resolution element is assumed to 
be proportional to the number of stars with line-of-sight ve- 
locities corresponding to this wavelength interval. The line- 
of-sight velocity distribution measured from the line profiles 
is a discretized function linearly related to the underlying 
DF [cf. eq. (1)]. This linearity is lost when line profile mea- 
surements are represented by the quantities v, a, h^, hi — 
these quantities are obtained by a least-square fit of a Gauss- 
Hermite series to the observed line profile. 

To re-express the observed kinematics in terms of quan- 
tities that depend linearly on / we proceed as follows. Con- 
sistent with the assumption of spherical symmetry, we set 
the mean streaming velocity v and all odd velocity profile 
moments to zero. Next, we obtain an estimate for the veloc- 
ity dispersion (second moment) a , by integrating over the 
line profile specified by (a, ^4); for negative /i4, until it first 
becomes zero. For small /14, the linear correction formula 
a = a{l + y/Qhi) holds ( van der Marel & Franx 1993), this 
correction results ina>a for peaked profiles with /14 > 0. 
The numerical correction from integrating over the velocity 
profile also has this property (BSG). 

From the measured hi{Ri), we compute new even 
Gauss-Hermite moments s„{Ri;a) by expanding the series 



£(7?„U||) ^CoY^ h,{R{)H,(x) exp(-a;V2) 



j=0,4 



{ho = 1, h\=h2 = h-j, = 0) with x = Vn/a{Ri) as 
C{Ri, V 11) ^ ^ A^„ s„ ii'„(2-) exp(--E^/2), 



(24) 



(25) 



n=0,2,4,... 



where x = v\\/a(Ri). Here Hn are Hermite polynomials, 
the Nn are normalization constants (G93), and a{Ri) are 
fiducial scaling velocities generally different from a{Ri). In 
practice, we have found it convenient to take for a{Ri) the 
velocity dispersions aiso{Ri) of the isotropic model in the 
given potential $(r) with the same stellar density as the 
galaxy being analyzed. The Sn{Ri;a) are estimates for the 
Gauss-Hermite moments related to the true velocity profile 
by 



,{R^;a) = (2"-in!) 



1 r\-l/2 



(26) 



For a theoretical model, we obtain corresponding mo- 
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ments by inserting the projected DF N{Ri,v\\) from (1) into 
eq. (26) instead of C[Ri,v{), using the same a(Ri). Clearly, 
the new Sn-moments of the composite model are linear in 
the ak, i.e., 




(27) 



with the 



Sn corresponding to the respective fk ■ 
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